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We show in a realistic d x 2_ y 2 symmetry gap model for a cuprate superconductor that the clean 
QQ ' vortex lattice has discontinuous structural transitions (at and near T=0), as a function of the 

magnetic field B along the c-axis. The transitions arise from the singular nonlocal and anisotropic 
susceptibility of the d x 2_ y 2 superconductor to the perturbation caused by supercurrents associated 
with vortices. The susceptibility, due to virtual Dirac quasiparticle-hole excitation, is calculated 
^ . carefully, and leads to a ground state transition for the triangular lattice from an orientation along 

I 1 one of the crystal axis to one at 45° to them, i.e, along the gap zero direction. The field scale is seen 

to be 5 Tesla ~ (Ao/£a) 2 $o, where Ao is the gap maximum, t is the nearest neighbour hopping, 
a is the lattice constant, and $o is the flux quantum. At much higher fields (~ 28T) there is a 
ryj ' discontinuous transition to a centred square structure. The source of the differences from existing 

, calculations, and experimental observability are discussed, the latter especially in view of the very 

small (a few degrees K per vortex) differences in the ground state energy. 
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An external magnetic field enters a (type II) superconductor as a collection of quantized magnetic flux tubes. The 
flux tubes with associated supercurrent vortices form a triangular lattice jjj as is seen in conventional superconductors 
' Deviations from this structure are of considerable interest. In conventional superconductors, the observed 

, deviations have been attributed to the anisotropy of the underlying one electron energy spectrum. In heavy fermion 
00 ' an d high T c superconductors, an additional and potentially very interesting reason is the existence of an unconventional 
0^ \ superconducting order parameter, with gaps which have nodes and change sign. Indeed there is experimental evidence 
both in heavy fermion systems [Ej the exotic superconductor SrRuC>4 jsj and in cuprate superconductors (|J^]that 
the vortex lattice is not triangular (for some field and temperature regimes). Recently, measurements of small angle 
neutron diffraction from untwinned YBa2Cu307_a single crystals shows a very well formed triangular lattice that 
undergoes an orientational transition from along a axis to along b axis for a 3T magnetic field at 33° to the c-axis 
||. The reasons for possible nontriangular structure as well as for the structural transition arc not fully established, 
and are the subject of considerable theoretical work P |l3]|. In high T c superconductors, entropic effects, abetted by 
high transition temperature as well as weak interlayer coupling play an important role, and the classical statistical 
I ■ mechanics of interacting, meandering flux lines, of the vortex fluid phase, and the solid fluid transition has developed 
' O I into a major theoretical and experimental subfield fl4|| . 

In this paper, we focus on the nature of the flux lattice at temperatures well below the superconducting transi- 
tion, where vortex configurational entropy effects mentioned above are negligible. The ground state structures and 
structural transition then directly reflect the electronic peculiarities of the superconductor, and thus probe the latter. 
For cuprate superconductors, a number of measurements show that the superconducting gap At has nodes |fl5|| , has 
a magnitude with a (Ao/2)| cosfc x a — cosfcj,a| dependence on the two dimensional wave vector k across the Fermi 
surface and that transport lifetimes of quasiparticles are long [|l7| for T <C T c . Thus one can assume well defined 
low energy, nodal quasiparticles, with an experimentally determined one electron dispersion fis) ] and gap function 
Aic fl^l - The question of interest is the effect of the zero gap, anisotropic Dirac like linear quasiparticle excitation 
spectrum on the interaction between vortices, and thus on vortex lattice structure. There is considerable evidence e.g., 
from magnetic field dependent electronic specific heat |l9| , electronic thermal conductivity |2(]] , and superfluid density 
fHf measurements that an external magnetic field going in as vortices has a strong effect on electronic states, changing 
their density and lifetime. The relevant issue here is somewhat the reverse, namely the effect of the quasiparticles on 
interaction between vortices. The order parameter phase associated with the vortex, and the related magnetic vector 
potential together constitute the superfluid velocity field v s (r)[= J^i v s(r — R-z)] where the vortices are located at 
points R;. The extra superfluid kinetic energy, being quadratic in v s , clearly has a part that depends on two vortex 
coordinates and is thus structure sensitive. In addition to this 'diamagnetic' term, which is the additional kinetic 
energy of the rigidly moving superfluid and which is minimized for a triangular lattice [1], there is another 'paramag- 
netic' term due to the perturbation of quasiparticles by the superfluid current via the term k • v s . This causes virtual 
particle hole excitations; two vortices interact via the exchange of quasiparticle quasihole pairs. This polarization 
term depends on the quasiparticle excitation spectrum. In a clean s-wave superconductor, the process leads to an 
additional isotropic interaction between vortices of order (H/H C2 ) relative to the diamagnetic term. However, in a 



d-wave superconductor where the excitation gap vanishes along some (nodal) directions, one expects the nonlocal 
polarizeability to be larger as well as anisotropic; this gives rise to an interaction between vortices which depends on 
the orientation of the line joining them with the respect to crystaline axes, and consequently can be the cause of novel 
long range positional order. 

The ground state energy arising from quasiparticle / hole mediated interaction between vortices depends linearly 
on the nonlocal current susceptibility X^(q), for wave vectors q equal to the reciprocal lattice vectors G of the 
vortex lattice. Because the gap as well as density of quasiparticle states vanish linearly near the node, X^fl(q) is 
proportional to \q x \ or \q y \ for small q. This nonanalytic behaviour, noticed first by Kosztin and Leggett [ p2[ , has also 
been discussed by Franz et al who were the first to analyze microscopically its effect, as well as of the anisotropy 
in Xa/3' on vortex lattice structure at T = 0. These authors found a rich phase diagram in the field temperature 
plane, with a centered rectangular lattice at T = whose inner angle varies continuously as a function of field, as 
well as a sudden orientational transition at higher temperatures, and a transition to a centered square lattice for very 
high fields and low T. In obtaining these results, Franz et al made a "local" approximation for the gap function, i.e., 
assumed Ak = Ak+c and more importantly they assumed a momentum independent quasiparticle current which 
leads to a response function Xxx(^) ~ The anisotropy then enters only through X™(G). We carry out here 

a more detailed and realistic calculation of the nonlocal susceptibility, considering the strong k dependence of the gap 
function Ak and quasiparticle current jk properly, and using a realistic one electron dispersion. The diagonal terms 
Xxx an( i Xyy are unequal and large, and this anisotropy is seen to be the underlying cause of the transition. The 
contibution of the off diagonal susceptitiblity Xxy ^ s smaller than that of the diagonal susceptibilities. Our results for 
structural stability (at T = Q) are therefore quite different from those of Franz et al 0,^3) ■ 

Confining ourselves to T = 0, we find, as summarized in a phase diagram (see Fig. ^7 that the stable lattice at low 
fields is triangular. At about 5 Tesla (for the parameters chosen) the orientation of the smallest G vector changes 
from being along one of the axes to being along the order parameter node direction, because the system is most 
susceptible to excitations with wave vector along the node. We have analyzed the driving force for this transition, 
both analytically and numerically, and find that it arises from a subtle balance between the term linear in |G| and the 
quadratic term, which are slightly different for the two orientations. The field scale for the transition is approximately 
given by the condition (taG/ Aq) ~ 1 which is natural on dimensional grounds. The Fermi velocity is ta and the 
energy scale associated with the superfiow (V0) with Fourier component G is therefore taG. The polarizability 
or susceptibility has an energy scale (1/Ao) where Ao is the gap magnitude which sets the scale for quasiparticle 
excitation energies. Thus the dimensionless susceptibility of interest is (taG/ Aq). For realistic parameters t, a, and 
A this translates [(taG c / A ) ~ 0.37] to a field scale of 5.2T. 

We find that the node oriented triangular lattice is stable till about 28T, whereupon a discontinuous transition to a 
centered square lattice takes place. This structure, which is orientationally commensurate with the symmetry of the 
quasiparticle dispersion, is probably the most stable T = phase when electronic commensurate effects dominate. 
However, the calculated field scale is large enough that the London approximation used, valid for H <C H C2 , is not 
reliable, because vortex core effects cannot be neglected at these high fields. 

In the next section (Section II), we describe the model and the theoretical approach used. The tight binding 
quasiparticle Hamiltonian is decomposed into an unperturbed part Hq and a term Hj due to the quasiparticle vortex 
interaction. The free energy or the ground state energy can be obtained as a power series in Hi or equivalently the 
density of vortices. For low vortex densities (H <C H C2 ) the leading or n\ term is sufficient, and describes quasiparticle- 
hole mediated vortex interaction, in addition to the bare superfluid kinetic energy. We discuss the former carefully 
in terms of the nonlocal, anisotropic current susceptibility XaflCl) smce t ne energy can be expressed as a reciprocal 
lattice vector sum over X„ /3 (G). We obtain x^o(q) semianalytically for small q at T = 0, as well as numerically 
(Section III). 

The calculations for different two dimensional structures are discussed in section III. For a given magnetic field B 
the most general centered rectangular lattice [01,03] can be described in terms of an angle 9 related to the aspect 
ratio (01/02) as tan# = 01/02, and an orientation <p with respect to crystal axes (Fig. |^). We compute the ground 
state energy as a function of these two variables for different magnetic fields. The basic vortex related electronic 
energy parameters are the following. The vortex has three energy scales, namely the diamagnetic single vortex 
energy, of order 3450K per vortex, the vortex interaction energy, the diamagnetic part of which has a value ~ 1440.fr 
(for nearest neighbor vortices), and the paramagnetic interaction term which is about 3bOK at 5T field. The last, 
and the smallest term is structure sensitive and is of interest here. The ground state is analyzed as a function of 
8, 4> for several field values in section III. It turns out that the structure sensitive part of the last (paramagnetic 
vortex interaction) term is extremely small, of order a few degrees per vortex. This has obvious implications for the 
observability of the transition, because the structural changes predicted and the clean limit anisotropies obtained can 
be easily overwhelmed by effects of disorder, eg., vortex pinning and the muting of the paramagnetic susceptibility 
anisotropy, and nonanalyticity by disorder. However, the size of the structure sensitive terms is larger, the greater 



the (vf/va) ratio or anisotropy. One can thus imagine situations where this effect is quite iarge. 

In the concluding section (Section IV) we briefly discuss thermal effects, the consequence of the predicted transition 
and their observability, the reason for which our result doffers from earlier, and the experimentally observed structural 
transitions. 



II. THEORY 



A. Model 



We consider a two dimensional lattice model with nearest neighbor and next nearest neighbor hopping for a Cu02 
plane of high T c superconductors. The (mean field) Hamiltonian in this model is given by 

H = -t ( c \<j c j° + h - c -) +*' + h - c -) + ( A ^ c lT c ]i + h - c -) - vJ2 c ™ c ™ ' w 

<ij>\o <<ij»io <ij> i 

where t and t' are nearest neighbor and next nearest neighbor hopping integrals respectively. This corresponds for 
appropriate choices oft and t' to an open Fermi surface which is observed in angle resolved photo emission experiments. 
The pair amplitude Ay is considered to be d x 2_ y 2-wa\e like, i.e., A.i t i± a x = — Ai^± a y, where a is the lattice constant 
in a square lattice, and /U is the chemical potential. 

When we apply a magnetic field beyond the lower critical field H Cl in high T c superconductors, the magnetic field 
goes into the system in the form of vortices. The magnetic induction is screened over a length A, the penetration 
depth. The pair amplitude acquires a phase: Ay — ► Ay exp[— where Oij is the sum of polar angles of all the 
vortices measured with respect to a particular axis, for the centre of mass of the pair ij. We write y as (9i + 9j)/2, 
(as an average of the angles of individual coordinates of the Cooper pairs), which is consistent upto 0(l/kp{;) 2 , where 
hp is the Fermi momentum and £ is the superconductive coherence length. The pair amplitude (order parameter 
magnitude) vanishes at the center of the core of a vortex, and over a distance £ it acquires its uniform value. For a 
collection of vortices with H <C H C2 , i.e., with intervortex spacing 3> £, or the London limit, we assume the order 
parameter magnitude be uniform, throughout the superconductor (there are 5 function sources of phase rotation at 
the locations of vortices). There is a vector potential A such that V x A(r) = B(r) where B(r) is the local magnetic 
induction, along the c-axis. Its effect in this model is to change the hopping integrals to 

(t, t') -» (t, t') exp[i(e/hc) [ J A ■ dl] (2) 



for hopping from site j to site i. We then make a gauge transformation: c, CT — > Ci a e l0i / 2 . We thus obtain the 
Hamiltonian as 

H = -t Y ( C L C J° ex P[*(^ -0j)/2 + i(e/hc) / A • dl] + h.c.) + ^ ( A v4] c ]l + h - c ^ 

<ij>io r * <ij>ia 

+t' E (4 c *r "WW* - °j)l 2 + *( e A c ) / ' A • rfl] + h.c.) C *W • (3) 

<<ij>> ri i 

The phase difference between two nearest or next nearest neighbour sites can be expressed as 



- 6j) + (e/hc) J 3 A ■ dl ~ (n - Tj) ■ (U/ t 9 ~ (e/fic)A.) = (m/»)(r< - r,-) • v s (r 2 ) . (4) 
Here the superfluid velocity 



v s (r) = — 
m 



(5) 



for a single vortex and for a collection of vortices, v s (r) = J^i v s(r — R-z) (where the vortices are located at R/. We 
then assume that the phase difference between two neighboring lattice sites is very small (which is certainly true in 
the London limit) so that we expand exponentials in Eq. (^) upto quadratic terms. 

Using Eqs. (5) and (4) in Eq.(3) for H, and expanding upto quadratic order in the small quantity v s (r) we have 



H = Ho+H I + Hu, (6) 

where the free Hamiltonian 

Ho = J2 & C L C ^ + H[ A k c k T c -ki + h - c -} . ( 7 ) 

k,cr k 

with £k = — 2t[cos(k x a) + cos(k y a)] + At' cos(k x a) cos(k y a) — fx and A k = (A /2)[cos (k x a) — cos (k y a)], A being the 
maximum quasiparticle excitation gap. Here k lies in the first atomic Briliouin zone, i.e., — ^ (k X: ky) — . A 
typical structure of the Fermi surface is shown in Fig. ||. Gapless quasiparticle excitations exist along k x = ±k y 
directions as noted in the figure. The interaction term (first order in v s ) can now be expressed as 



c k(T c k+G(T - c k+Go .c ko -] [mv*(G) sm(k x a) + mv^(G) sin(k y a)} 

k G>0 

-4(at'/K) ^ E [ c L c k+Gcr - 4+Ga c k"] [ m v s ( G ) sin(k x a) cos(k y a) 



k G>0 

+mt;g(G)siii(fe v a) cos(k x a)) 

E EE VkG \- C la C *+Ga - cl +Ga C ka ] (8) 
kcr G>0 

Here Vk,G is purely imaginary. The term Hn is quadratic in v B and contributes to the free energy as a diamagnetic 
term. It is given by 

H n = 2N° (t - t') (a 2 /^ 2 ) ^ m <(G) m<(-G) (9) 

G 

with N® being the number of superfluid carriers, and a refers to cartesian variables x and y. Paramekanti et al. 
p|| have recently shown that the quantum phase fluctuation of the order parameter reduces the superfluid density 
considerably. We thus reexpress the term Hn phenomcnologically in terms of the measured A as 

V ' G 

where d is the mean interlayer separation of weakly coupled superconducting layers and A is the area of the system. 



B. Free Energy 

We now calculate the free energy as a power series in v s (r) or equivalently the vortex density. The diamagnetic 
or Ginzburg-Landau term, of first order in Hn, is the largest, and the structure sensitive part of it is known to be 
minimized for a triangular lattice (1). The energy does not depend on its orientation with respect to the crystal 
lattice. We are interested here additionally in the paramagnetic term, of second order in Hn Including this, and the 
magnetic field energy contribution, the free energy to second order in vortex density is given (per unit length along 
the c axis) by 



2Adh z c 



where G is the reciprocal vector of the vortex lattice. The individual vortex energy is not included here, as it is 
not relevant for the question of vortex lattice structure. x d — (c 2 ?i 2 d/47re 2 A 2 ) is the diamagnetic term arising from 
the term Hji (Eq.|lO|) to first order, and Xaa(G) is the paramagnetic current susceptibility due to second order 
contribution from Hj. Higher order contributions are neglected since the expansion parameter is (n v /n) where n v is 
the vortex density and n is the electron density. This ratio is obviously much smaller than one. 
The paramagnetic susceptibility xLg(q) is expressed as 



1 



(2tt)3 



tz / a pTT / a 



7r / a J —-R I a 



xlM = 17Ty i \ dk x &, 7Q/3 (k)n(k,q). (12) 



The current operator dependent terms 7 Q ^(k) are explicitly given as 



7xx(k) = [2as'm(k x a)(t - 2t' cos(k y a))] , 
7 yy (k) = [2asm(k y a)(t — 2t' cos(k x a))] 2 , 

7xy(k) = [2a sm(k x a)(t — 2t' cos(k y a))] [2asm(k y a)(t — 2t' cos(k x a))] 
= 72/* (k) • 



(13a) 
(13b) 
(13c) 
(13d) 



The zero frequency susceptibility of wave vector q for quasiparticle quasihole of momentum k is II(k, q) and has the 
form 



n(k, q ) = 



1 



1 



CkCk+q + A k A k+q 



E^E^ 



(14) 



with the quasiparticle energy = 



It is expected that the above susceptibilities are anisotropic due to the nonlocal nature of Ay, reflected in the k 
dependence of A k . Though anisotropic, they possess certain symmetries: Xaailx, ~%) — Xaailx, %) = Xaa(~1x, Qy)', 
Xx V (qx,-q y ) = -Xxy(qx,q v ) = Xly(-Qx,q y ), and x xx (qx,q v ) = Xyy(Qy,Qx)- These symmetries suggest that the 
susceptibilities are functions of \q x \, \q y \ and sign{q x q y ) only. A naive perturbative expansion of X^«(q) in powers of 
q fails since the coefficient of quadratic term in q is divergent, due to the vanishing of A k on the Fermi surface at 
k x = ±k y points. We however proceed to evaluate these analytically as follows. 

We write 



Xlp(qx,q y ) =^2Xap(Qx,qy) 



(15) 



where XaR ^ s ^ ne contripution of j-th quadrant (j =1-4) of fc-space. For instance, 



(27T) 



7T / a 



ir/a 



dk y -f a(3 (k)U(k, q) 



(16) 



is the contribution due to 1st quadrant. We present the calculation of X^aiQx, Qy) below in detail. 

In terms of an alternative coordinate system (ki, k%) whose origin is at the nodal point on the Fermi surface as shown 
in Fig. |J, the old coordinates in the first quadrant are expressed as k x — -^(ko + k\ — k^) and k y = -^(ko + k\ + k%), 

where ko is defined as fi = — At cos(fc a/v / 2) + At' cos 2 (fc a/v / 2). We use £ k ~ fovpki and A k w %VAk% in linear form, 
where the Fermi velocity vf = (^75)^ _ 2f' cos sin and va = (^75)^0 sin ^| is the velocity of quasiparticles 
along the &2 direction. Since in ci-wave superconductors vp S> va, the phase space of ki effectively is much larger 
than that of k\ for a given value of quasiparticle energy. We observe that this is the cause of strong anisotropy in 
the diagonal susceptibilities as we see below. If <j> be the angle of a k-vector with k% axis in this new coordinate 
system, £ k ~ i? k cos</> and A k ~ i? k sin^>. By Taylor expansion in q we find £ k £ k + q + A k A k+q ~ E^E^ + a q ) and 
E^Ek+q ~ \E^ + E'lcCkq + /3 q |, where a q = H(qiVF cos0 + g2«A sin</>) and /3 q = (ft 2 /2){q\Vp sin^ — qiVA cos<?!>) 2 with 
91,2 = {q y ± q x )/V2 respectively. Since a q is negative for some region of (f>, the quantity (i? k + E^a^ + /3 q ) may 
be negative as well as positive which we refer below as the region I and II respectively. It is negative in the regime 
Ei < E k < E%, where E® w (— /3 q /a q ) and E% « (— a q + /3 2 /a q ). Expanding 7 a /s(k) upto linear in £ k , we perform 
the integrals over i? k in for both regions I and II separately to obtain 



Xxx (qx,q y ) 



VfVa 



Xy y (qx,q y ) 



ir 2 h 2 VFVA 



J # {fti (-a, + 2(/3 2 /a q ) In |a q /A |) 
+ ^(2?+ ft 2 +2^ft 3 ) (a 2 - 3/3 2 - 3/3 2 In |a q /A„|) 
+ / <*0 { 21n ( 2 ) (/3q/«q)^i - 2/3 2 (P+^ 2 +P^^ 3 )ln(a q /A )} 
y # {7l x (-a q + 2(/3 2 /a q ) In |a q /A |) 
+- |(X»J ft 3 + T>+ n 2 ) (a 2 - 3/3 2 - 3/3 2 In |a q /A |) | 



(17a) 



Xxy(9x,q y ) 



^ i<f> {21n(2) (fc/aJHi - 2f3l(V+K 3 + Z£ ft 2 )ln(a q /A )} 
^ rf0 {^i (-a q + In |a q /A |) 

&R2 + Ka) cos (a q - 3/3 q - 3/3 2 In \aj A 1) 



(17b) 



17 



{21n(2) {ft/ajKt - 2/? q (^ 2 + ^ 3 ) cos ln(a q /A )} 



where 2?J = cos ^ ± (vf/va) sin < 



= (i 2 + /it') sin 2 (fc a/\/2) , 
K 2 = 2V2(-^—)t'(t - 2t'cos(/coa/v / 2))sin 3 (fc a/\/2) , 

TlVF 

TZ 3 = V2(-^—)(t 2 + fit') sm(k a/V2) cos(k a/V2) , 
Tivf 



(17c) 



(18a) 
(18b) 

(18c) 



and Jj and Jjj represent the integrals over 4> < 4> < 2ir) for the regime of <f> in which a q < and > respectively. 
We, similarly, calculate y^i for three other quadrants. We observe that Xxx and \yy differ substantially through the 
last terms within both the integrals Jj and J J7 in their expressions (17a) and ( |l7b| ) since vp ^> fA- These lead to 
an anisotropic diagonal susceptibility. We note that these terms also the terms involving J„ arise due to keeping the 
linear dependences of k\ and fc 2 in 7 Q| g(k). However for an approximation k x = k y = fco/v2 in j a p(k), Xxx — Xyy as 
obtained by Franz et al [Q . 

Since the angular integrals in the expressions for X^g cannot be performed analytically, we numerically integrate 
these to obtain x V a g in the next section. We shall then compare these semianalytically obtained X V a g with those 
completely numerically obtained through Eqs.(jl^)-(|l4l). 

We now turn to obtain the equation for the vector potential in the gauge G Ag = (from Eq. jll|)). By minimizing 
the energy with respect to as Ac, we have 



(A G ) a 



Aire 
TicG 2 d 



[ x d 5 aB -xl a (G)\ 



We thus obtain 



where 



Qa^(q) = jjS aP - [ ) x%{<\) 



Aire 2 



Using Eq. 



in Eq. ( 1 lh , we get 



8irA 



1 



c 2 h z d 



GaQafiGg 



DetQ 



B. 



G 



We now express B G in terms of N v ,$>0i Qap and G. For a vortex lattice, 

(V0) G =2mN v G ^ ' 



G 2 ' 

where N v is the total number of vortices. We then obtain from Eqs. (|l^) and ( f20| 

Det Q + Qape^e^G^Gs 



B G = N v $ 



G 4 + G 2 Q aa 6 a3 + DctQ 



(19) 



(20) 



(21) 



(22) 



(23) 



(24) 



with e 12 = 1 = — e 21 , e 11 = = e 22 , and <&o = hc/2e is the quantum of flux. Therefore, the free energy for a vortex 
lattice per unit volume becomes 



G 



Det Q + Q a f3e aj e l3S GyGs 



G 4 + G 2 Q aP 5 aP + Det Q 



GaQafiGf: 

DctQ 



This has an approximate but much simpler form as 



1 (X> \2 G^Qyy 



8?r 



G 



G 4 



(25) 



(26) 



which is essentially important for determining the ground state structure of the vortex lattice. This form is exact when 
|§(V0)g| 3> |eAc|/c which is true. The component G = will give free energy T for average magnetic induction B. 
For determination of vortex lattice structure, one should in principle minimize Gibbs free energy Q = J- — BH/An. 
Here H is the applied magnetic field. Beyond H Cl , magnetic field penetrates the superconductor almost fully. Thus 
B ~ H, specially so in high T c superconductors, since H Cl <C H C2 . B does not vary much for different vortex lattice 
structures for a given H as we see in our numerical study that the ratio (H — B)/B ~ 1CP 7 . We, therefore, minimize 
AJ 7 = T — T ie. that part of the free energy which depends on G, in Eq. ( j26| ) for different choices of nonzero G's 
corresponding to different structures and with a cutoff G < 7r/£. 



III. NUMERICAL STUDY 



The values of the phenomenological parameters that we have used for the numerical computation of Q a f3(Qx,Qy) 
and later for the free energy are taken from angle resolved photoemission experiments [[l6 18 , penetration depth 
measurement and band structure calculations J25J for high-T c compounds. These are as follows: t = 1150K, 
t' = 0.48*, A = 400K, a = 3.8l, d = 10A, £ = 20A, A = 16001, and (j, = -1.33* which corresponds to doping 
concentration x ~ 0.19. 

Using standard Gaussian quadrature, we integrate over k x and k y in Eq. (|l^) to obtain X^ (3 (q). In Fig. |s| we have 

shown the dependence of paramagnetic susceptibilities (a) Xxx > O 3 ) Xy V > an d (c) \xy m units of \ d f° r positive q x at 
different positive values of q y . Susceptibilities for negative values of q x and q y can be obtained by using the symmetries 
discussed in the previous section. We see that Xxx(lx,qy) 7^ Xyy(<lx,£[y) in general. This strong anisotropy i n th e 
diagonal susceptibilities is due to the strong k dependence of nature of Ak and the k-dependent 7 Q/ g(k) (13a)-(13d). 
The diagonal susceptibilities are large compared to off-diagonal one 



We numerically fit, guided by the semianalytical form in Eqs.(17a)-(17c), to obtain the approximate functional 



form of XapiQx* 1y) for fea, q y a < 0.3 as 



(Qx,q v ) 



> 



%\ 



0.31(5|flf x |a) + O.U(Sqxa) 2 - 0.35(5 fe a) 2 In \Sq x a\] if \q x 

0.35(S\q y \a) - O.U(6q y a) 2 + O.lO^a) 2 In \Sq y a\ 

(0.10 + ^- a ){8q x af + (-0.16 + ^L.)(Sq x a) 2 In \Sq x a\\ if \q x \ < \q y 



X yy (Qx, qy) — X?a;(9j/j 1x) i 



Xl y {qx,q y ) =7 



(0.11 + —) q< + (0.15 



q> 



0.18. 

q> 



q\ + (0.07- 



0.12. 

q> 



sign(q x q y ) 



(27a) 

(27b) 
(27c) 



with g>,< = max, mindg^l, | qr y | 7 = ^f( ^§§? )t whose numerical value is 1.18, and the parameter S — t/Ao. These 
phenomenological forms can be explained from the semianalytical expressions ( p_7a) )-(17c) as follows, (i) Firstly, why 
do Xxxilx, q y ) and Xy y {qx,q v ) not depend on the signs of q x and q y , and Xxy(qx7qy) does depend on sign(q x q y )? This 
is due to the symmetry reason discussed following Eq.(pT[). (ii) Why does x xx (qx, qy) does mainly d epen d on whether 
\Qx\ > Iffy I ° r not? This can be understood by the following exercise. We find a term from Eq.( 17a ) as \q y + q x \ 
assuming q x ,q y > 0. The corresponding term will be \q y — q x \ when we consider the contribution from 2nd quadrant. 
When we add these two contributions, we see that the sum depends on the greater of q x and q y . (iii) Following the 
argument above in (ii), the difference between the two terms is smaller of q x and q y . This is the reason why Xxyfax, qy) 
depends mainly on the smaller of \q x \ and \q y \. (iv), Since Xxx(qx,q y ) ^ Xyy(Qx,Qy), and xL(fe,%) = xl y {q y ,qx) 
for symmetry reasons, the dependence of X xx (qx-, qy) on q x a nd q y is asy mmetric, (v) The linear, quadratic, and the 
logarithmic dependences on q follow from in the expressions (17a )— ( |17q ) . 



We next numerically perform angular integrals in equations (17a)-(18c) along with the contributions from other 
three quadrants to obtain semianalytical susceptibilities and then compare with the fully numerically obtained sus- 
ceptibilities. In Fig. ^| we have shown Xxx(lx,0) and XyyiQx,®) evaluated in the two ways. The linear approximation 
of energies in the analytical expressions is a good approximation for determining linear dependence on q x as we see 
in Fig. that they agree for very low q x . They however differ for higher q x since our analytical expressions are not 
consistent in determining quadratic dependences on q as we have neglected higher order k dependences to quasiparticle 
energy. It is however clear that x xx 7^ Xyy which is our main result. 

We consider a face centered rectangular vortex lattice (as shown in Fig. ||) with area of the unit cell A = 2$o/F, in 
general. Angle 9 determines the sides of the rectangle with a fixed area. The sides of the rectangle are ai = [Atanfl] 1 / 2 
and a 2 = [A/ tan^] 1 / 2 . We then readily obtain reciprocal lattice vectors for a vortex lattice, in general, to be 

G mn (B,9) = (n + m) — e x + (n - to) — e y , (28) 
oi a 2 

where n and to are integers (both positive and negative) including zero. If the vortex lattice makes an angle </> with 
the underlying atomic lattice, we find 



G m n(B, 9, 



+e y 



( n + to) — cos — [n — to) — sin c 
a\ a 2 

, s2tt . , . ,2tt 
(n + to) — sin <p + (n — to) — cos ( 
ax a 2 



(29) 



The lattice is a centered square for 9 = 45° and triangular when 9 = 60°. There is symmetry of rotation about 
= 45°, since the lattice is considered as centered rectangular. We therefore need to determine free energy for 
45° < 9 < 90° and < < 45°. 

We then numerically compute the free energy per vortex (without the single vortex energy which does not depend on 
structure), AF — dAT jn v using Eqs.(p6[) and ( |29| ) as functions of the parameters 9, 1 , and B. The reciprocal lattice 
vector G changes with the change of any one or more of the parameters. Thus AF differs for different combination 
of these parameters (<fi, 9, H). In Fig. [t] we have shown the dependence of AF at a low field B = 2 Tesla as a function 
of <j) for the angles 9 — 60° and two neighboring angles 9 — 58° and 62° (on either side of 9 = 60°). It is clear that 
AF is a minimum for the triangular lattice. We notice also that AF is a minimum for the triangular lattice when 
<f> = 0° and 30° which in fact correspond to the same lattice configuration. We likewise find that in the whole of the 
low field regime, the ground state configuration of the vortex lattice is triangular with one of its arms parallel to one 
of the crystal axes. 

Interestingly, the orientation of the lattice changes discontinuously as we increase the magnetic field though the 
structure continues to be triangular. In Fig. || we have shown the dependence of AF on cj> for the triangular lattice 
configuration at three chosen fields 2, 5, and 8 Tesla. At nearly about 5 Tesla field, AF is minimum for all 0°, 
30°, 15°, and 45° orientations; the latter two angles correspond to same lattice configuration, like the former two 
angles. On the other hand at the field of 8 Tesla, AF is minimum for <p = 15° and 45° only. The triangular vortex 
lattice changes its orientation discontinuously at about 5 Tesla field. While the triangular lattice has one of its 
arms parallel to one of the crystal axes at lower field, it aligns to one of the crystal axes by 45° at higher field. 
We understand this discontinuous transition by comparing the energies contributed to AF by the G vectors of the 
lowest magnitude (since they contribute most to the free energy) for these two preferred orientations. Considering the 
symmetries of the susceptibilities, it is sufficient that we consider only those G vectors which have positive G x . We 
thus consider 3 G vectors for each of these two orientations. These are (a) (1/2,±V3/2)G and (1,0)G for = 0° and 
(b) (V3- 1, V3 + 1)G, ^(V3+ 1, V3- 1)G, and (1/V2, 1/V2)G for = 45°, where the length of the smallest G 

vector G = 27r(2/ v / 3) 1/2 (F/$ ) 1/2 - In Fig. g we have shown the energy E contributed by these individual G vectors 
to AF for 2, 5, and 8 Tesla fields. We find the total energy contributed by above 3 G vectors for <fi = 0° and 45° 
orientations as (a) 296.05K and 296.51K for B = 2 Tesla, (b) 282.96K and 283.11K for B = 5 Tesla, and (c) 274.15K 
and 273. 60K for B — 8 Tesla respectively. Clearly, the triangular lattice makes an orientational transition at about 5 
Tesla field. 

To understand the field scale 5 Tesla for the above orientational transition, we compare the energy contributed by 
the above 3 G vectors for each of the preferred orientations. The ratio of these energies can be expressed as a function 
of a = tGa/Ao using equations d26|)-(|29|). This is given by 

F 2 3 -/ 2 (o;) 



where E\{E2) is the energy contributed by the above corresponding 3 G vectors of triangular lattice with 0°(45°) 
orientation. In Fig. [l0[ we have shown the ratio /i (a)/f2 (a) as a function of a. The orientational transition takes 
place when the ratio is unity. This corresponds to a c = tG c a/ Aq ~ 0.37. Therefore the critical field at which the 
transition take place, B 1 ~ (0.37/2\/27r) 2 \/3(A /at) 2 $o - 5.2T. 

The structure of the vortex lattice remains triangular with 45° orientation to the crystal lattice as we have shown 
AF in Fig. [fl] for a field as high as 25 Tesla. However, it makes yet another discontinuous structural transition to a 
centered square lattice with its axes parallel to the crystal axes at yet another critical field B 2 whose value is about 



28 tesla. Figure 12 shows that AF is minimum for 9 = 45° and <j> — 0° at B — 28 Tesla. The overall phase diagram 
for the ground state of the vortex lattice structure at T = is shown in figure [l| 

IV. CONCLUDING REMARKS 

We conclude by briefly discussing a number of questions such as the nature of the approximations used, the effect of 
nonzero temperature, consequences of the transitions, their observability, the reason why our conclusions differ from 
those found earlier, and the structural transitions experimentally observed. 

We have calculated the ground state energy assuming effectively that the interaction between two vortices is 
unaffected by the presence of other vortices. This is obviously a low vortex density approximation which seems quite 
reasonable since the dimensionless ratio (n v /n) is about 1/2500 for a field of a Tesla. However, we have not calculated 
the higher order corrections which while nominally of higher order in (n v /n) might have large or even divergent 
coefficients. Since the vortex vortex interaction depends on quasiparticle-quasihole susceptibility, a change in their 
spectrum due to the supercurrent (the Volovik effect p6| ) could have serious consequences. Here, we would like to 
make two points. Firstly, a calculation by Amin, Affleck and Franz [ fL3f , using a semiclassical approximation to include 
the nonlinear effect of the magnetic field a la Volovik, finds that this has little effect on the structural transformations 



calculated by them. Secondly, all the recent fully quantum mechanical calculations []27 28 1 of the density of Dirac 
quasiparticle states in a vortex lattice in the London limit find that for Bravais lattices, the density of states vanishes 
linearly with energy as in the absence of a magnetic field; only the slope changes. A general argument for this has 
been presented by Vishwanath p9fl . For these two reasons, we believe that our low density approximation is reliable. 

In the London approximation, the vortex cores are treated as S functions. In reality, they have a width of order 
the coherence length. We believe that the consequences of this approximation, at least for the low field structural 
transition, are small. The reason is that the structure sensitive part of the energy arises from the difference in 
contribution of the smallest reciprocal lattice vectors (Section III). For these, G£ ~ (1/10) at typical magnetic fields, 
so that the phenomenological assumption of a Gaussian vortex core with width ~ £ will make a negligible difference 
to the structure sensitive part. 

We have calculated only the ground state energy of the vortex lattice in this paper. At any nonzero temperature, 
there are obviously entropic contributions which could change the magnetic field at which the structural transition 
occurs, as a function of temperature. Here, we note that since both the structures (below and above 5T) are identical 
(triangular), and the structure difference sensitive part of the energy is a tiny fraction (< (1/1000) 2 ) of the vortex 
interaction energy, the elastic fluctuations in both structures are expected to be identical to order (I/I000) 2 so that the 
transition field should not be affected by temperature, so long as the input parameters (e.g., £k, Ak, A ) do not change 
with T. The same cannot be said of the high field (~ 28T) triangular to centred square lattice transition, because one 
has a tight packed structure and the other not. The expectation is that the former has lesser elastic fluctuations than 
the latter, so that the transition field boundary should shift to lower values with increasing temperature. However, 
this conclusion is tempered by the fact that the London approximation is unreliable at these high fields when vortex 
cores get close to each other so that our basic result may not be that reliable. 

One interesting consequence of the orientational transition at 5T, which might be measurable, is the change in the 
very low energy density of quasiparticle states. At least for a square lattice, Vishwanath |^9| has shown that there are 
quasiparticle states with linear dispersion and that there is a very small gap arising from higher order terms in the 
quasiparticle velocity. If this kind of result carries through for a triangular lattice, then it might be an experimental 
way of observing the transition. 

We have calculated here, for the first time, the actual energy of the structure sensitive term (Section III) and have 
found it to be small, of order a few degrees per vortex. Because of this reason, the transition might be difficult to 
observe, since pinning energies of larger size are generally present |3C|] , unless the system is extremely perfect. 

We have discussed in detail (Section II and III) the reason why our results differ from those obtained earlier. 
Basically, it has to do with the anisotropy of the nonlocal current current susceptibility, i.e. the fact that ~x^. x (G) ^ 
Xy y {G). The reason for this essentially is that we have an anisotropic superconductor. The \xx an< ^ Xyy functions 
are plotted in Fig.|[ 



The question of a non triangular structure of the vortex lattice in cuprates has attracted considerable experimental 
attention [^-^, especially since it has become established that they are d x 2_ y 2 superconductors. Earlier small angle 
neutron scattering measurements were on highly twinned 123 crystals, so that the observation of four fold diffraction 
symmetry does not imply rectanglar /square lattice. Moreover, the positional order is very poor. A more recent 
experiment Q on untwinned 123 single crystals shows much better translational order (higher G peaks are resolved 
for the first time) and a triangular lattice with axes oriented along a, distorted because of a-b asymmetry. The authors 
find no structural transitions upto 4T with field along c-axis. They however find a transition from a triangular lattice 
oriented along a to one oriented along the b axis at a field of about 3.8T, oriented at 33° to the c axis. This is certainly 
quite different from the transition to a triangular lattice at 45° to a axis at 5T predicted by us. As Johnson et al point 
out j^] the observed transition could be due to the presence of chains in 123, and the novel ab-c anisotropy caused by 
it (which may have a strong effect on many physical properties). In order to seriously explore our conclusions, one 
needs to do experiments on cuprates without chains and ideally with tetragonal symmetry, as again mentioned by 
Johnson et al ||] . 
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FIG. 1. Phase diagram for the structure of the vortex lattice. The position of the vortices are denoted by black filled circles. 
Bi and B2 are the fields at which the structural transitions take place as described in the text. The structures are shown 
diagrammatically. 





FIG. 3. A typical open Fermi surface for high-T c compounds which can be parametrized by a t — t' model. The shaded 
region denotes the occupied states with concentration 1 — x, where x is the doping concentration. The superconducting state 
is gapless along the diagonals of the Brillioun zone 
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FIG. 4. New coordinate system (ki, fe) in the first quadrant of the atomic BZ is shown. Its origin is at the 'nodal' point at 
which the superconducting gap vanishes on the Fermi surface. The length of the nodal vector is ko- 




FIG. 5. Dimensionless susceptibilities(a) \xx, (b) Xyy an d (c) Xxy ( m units of Xd) plotted against small positive q x a. The 
numbers associated with each curve are the corresponding values of q y a. The susceptibilities for negative values of q x and q y 
are related to the same for positive values of q x and q y by the symmetries discussed in the text. 
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FIG. 7. Free energy per vortex as a function of <j> for 9 = 58°, 60° and 62° at 2 Tesla field. 
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FIG. 8. Free energy per vortex for triangular lattice structure (8 = 60°) as a function of its orientation angles (j> for (a) 
B = 2, (b) B = 5, and (c) S = 8 Tesla fields. 
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FIG. 9. The contributions to energy per vortex by the three lowest G vectors of equal length for positive G x at three different 
fields for two different commensurate orientations of triangular lattice. The open(closed) symbols represents 0°(45°) orientation 
of vortex lattice with respect to crystalline lattice. The open symbols with a dot in their centers correspond to energy for two 
different G vectors with same G x . The circles, squares and diamonds correspond to 2, 5, and 8 Tesla fields respectively. 
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FIG. 10. The dashed line represents the ratio /i (a)//2 (a) as a function of a. The solid line is a guide to the eyes for the 
value of the ratio 1.0. 
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FIG. 11. Free energy per vortex as a function of 4> for 9 = 60° and 45° at a high field B — 25 Tesla. 
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FIG. 12. Free energy per vortex as a function of for 9 = 45°, 50°, 55°, and 60° at B = 28T field. Clearly AF is minimum 
for the structure corresponding th 9 = 45° and = 0°. 



